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Abstract 

We present a two-dimensional, well-balanced, central-upwind scheme for approximating 
solutions of the shallow" w-ater equations in the presence of a stationary bottom topography 
on triangular meshes. 

Our starting point is the recent central scheme of Kurganov and Petrova (KP) for 
approximating solutions of conservation law's on triangular meshes. In order to extend this 
scheme from systems of conservation laws to systems of balance law's one has to find an 
appropriate discretization of the source terms. We first show that for general triangulations 
there is no discretization of the source terms that corresponds to a w'ell-balanced form of 
the KP scheme. We then derive a new variant of a central scheme that can be balanced on 
triangular meshes. 

We note in passing that it is straightforward to extend the KP scheme to general un- 
structured conformal meshes. This extension allows us to recover our previous well-balanced 
scheme on Cartesian grids. We conclude with several simulations, verifying the second-order 
accuracy of our scheme as well as its well-balanced properties. 


Key words. Shallow water equations, central schemes, balance law's, unstructured grids. 
AMS(MOS) subject classification. Primary 65M06; secondary 76B15. 

1 Introduction 


We consider a flow in a. two-dimensional channel with a bottom elevation given by B(x), where 
x = (x, y). Let H(x, t ) represent the fluid depth above the bottom, and u(x, t ) = (u(x, t). v(x, t )) 
be the fluid velocity. The top surface at any time t is denoted by w(x, t) — B(x) + H(x. t ). 

The shallow water equations, introduced by Saint- Venant in [22], are commonly used to model 
flows in rivers or coastal areas. When written in terms of the top surface w and the momentum 
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(Hu, Hv ) these equations are of the form 
w t 4- (Hu) x + (Hv) y = 0, 
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( 1 . 1 ) 


This choice of variables is particularly suitable for dealing with stationary steady-state solutions 
(see [12, 21] for details). For simplicity we fix the gravitational constant, g, from now on to be 
9 = !• 

In this work we are interested in approximating solutions of (1.1) on triangular meshes. Our 
goal is to investigate how to adapt the semi-discrete central schemes on triangular meshes that 
were recently introduced by Kurganov and Petrova in [14] to this problem. We are interested 
in deriving a discretization of the source terms in (1.1) that preserves stationary steady-state 
solutions, as such solutions play an important role in the dynamics of (1.1). 

Central schemes for conservation laws have become popular in recent years as a tool for ap- 
proximating solutions for multi-dimensional systems of hyperbolic conservation laws. Like other 
Godunov- type schemes, central schemes are based on a three-step procedure: a reconstruction 
step in which an interpolant is reconstructed from previously computed cell-averages; an evolu- 
tion step in which this interpolant is evolved exactly in time according to the equations; and a 
projection step in which the solution is projected back to cell-averages. When compared with 
other methods, central schemes are particularly appealing since they do not require any Riemann 
solvers and systems can be solved component- wise. 

A first-order prototype of central schemes is the Lax-Friedrichs scheme [6]. A second-order 
extension is due to Nessyahu and Tadmor [19]. Extensions to two dimensions are due to Armin- 
jon, Jiang and Tadmor [1, 11], By estimating bounds on the local speeds of propagation of 
information from discontinuities, it is possible to pass to the semi-discrete limit (see [13, 15] and 
the references therein). There are several extensions of central schemes to unstructured grids. 
A fully discrete method is due to Arminjon et al. [2]; a recent semi-discrete scheme was pro- 
posed by Kurganov and Petrova in [14]. Balanced Central schemes for shallow water equations 
on Cartesian grids are due to Russo in the fully-discrete framework [21] (see also [18]) and to 
Kurganov and Levy in the semi-discrete framework [12]. 

There are many approaches to approximating solutions of (1.1). We refer, e.g. to [3, 4, 5, 7, 
9, 16, 17, 20] and the references therein. Our goal in this paper is to show that balancing is also 
possible with central schemes. We would like to emphasize that this is the first time' in which 
the balancing issues are treated for central schemes on unstructured grids. 

The paper is organized as follows. We start in Section 2 with a brief overview of the KP 
central scheme on triangular meshes. We note that this scheme is not limited to triangular 
meshes and it can be equally well applied to general unstructured grids. We also make the 
necessary adjustments to incorporate source terms into the scheme. We proceed in Section 3 
with the discretization of the cell-averages of the source terms for the shallow-water equations. 
The goal is to find a discretization such that the scheme will preserve stationary steady states, 
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i.e. zero velocities and a flat surface. We show that on general unstructured meshes, there is no 
discretization of the source terms in the shallow water equations that provide a well-balanced 
form of the KP scheme. We then proceed by showing how to modify the original scheme in 
such a way that it is possible to obtain a well-balanced discretization of the source terms. We 
conclude in Section 4 with numerical examples that demonstrate the accuracy of our scheme as 
well as its well-balanced property. 

Acknowledgment: The work of D. Levy was supported in part by the National Science Foun- 
dation under Career Grant No. D MS-0133511. 

2 Central Schemes for Balance Laws on Unstructured 
Grids 

W T e consider the two-dimensional balance law, 

u t + f(u) x + g(u) y = S(u.x,y,t), (2.1) 

subject to the initial data, u(x,y, 0) = uq{x. y). We are interested in approximating solutions 
of (2.1) that are computed in terms of cell averages on a fixed unstructured conformal grid. To 
simplify this exposition we first consider the scalar case. The scheme that is described below 
can be easily extended component- wise to systems of balance laws. We will make this obvious 
extension later on when dealing with the specific problem of the system of shallow' water equations 
( 1 . 1 ). 

We focus our attention on the central scheme on triangular meshes that was recently derived 
in [14]. We briefly overview the derivation of this scheme in the setup of conservation laws of the 
form 


u t + f (u) x + g(u) y — 0. (2-2) 

W r e assume a conformal triangulation of the domain consisting of cells 7) of area |7}|. The 
neighboring cells to 3} are denoted by Tj k , k = 1,2.3, while the edge that is joint between 7} 
and Tjk is denoted by E jk and is assumed to be of length h jk . We also denote the outward unit 
normal to Tj on the &-th edge as rijk, and denote the midpoint of E jk as Mj k (see Fig. 2.1). 

We assume that the cell averages on all the cells {Tj} are knowm at time t n , 

* j ±-\f T n(x,t n )dx, (2.3) 

and reconstruct a piecewise-polynomial 

u n (x,y) = y~'^u](x,y)x j (x : y). (2,4) 

J 

Here u”(x, y) is a two-dimensional polynomial that is yet to be determined and Xji x > y ) is the 
characteristic function of the cell Tj . To simplify' the notations we omit the time-dependence in 
uAx. v). We also denote hv n_-.Y or the Do'vnomial that is reconstructed in the cell T,> . 

j \ ~ J J a; ^ Xr J J 
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Figure 2.1: The triangular grid 


Discontinuities in the interpolant Uj along the edges of Tj propagate with a maximal inward 
velocity aj k and a maximal outward velocity a° v k . These velocities can be estimated (for convex 
fluxes) as 

af k (M jk ) = — min{Ai[J(uj(Mj*)) ■ n jk ], \i[J(u jk (M jk )) -n ifc ].0}, (2.5) 

a 'jk(M j k) = vaax{\ N [J{uj{M jk )) - n jk ], \ N [J {u jk (M j k )) -«jfc],0}, 

where J ( Uj ( Mj k )) is the Jacobian of the flux F — (/, g) evaluated at Mj k and < . . . < \n are 
its N eigenvalues. These local speeds of propagation are then used to determine evolution points 
that are away from the propagating discontinuities. An exact evolution of the reconstruction at 
these evolution points is followed by an intermediate piecewise polynomial reconstruction and 
finally projected back onto the original cells, providing the cell-averages at the next time-step 
n” +1 . Further details can be found in [14], 

A semi-discrete scheme is then obtained at the limit 


d_ n up'-u” 

— u ■ — lim — 


dt 


3 




At 


( 2 . 6 ) 


Most of the terms on the RHS of (2.6) vanish in the limit as At — > 0. The only quantity that 
has to be determined is a quadrature rule for the integrals of the flux functions over the edges 
of the cells. If we assume a Gaussian quadrature with m nodes 



m 

x)dx ss y^c s <p(x s ). 

5=1 


which is scaled to hj k . and denote the quadrature points on Ej k as G s jk , the KP scheme for 
triangular meshes is 


duj 

dt 


hj k 


IT 


aT n in . 4- 


fc=l 


a jk 


E [(a%F(u jk (G s jk ) + a$F( Uj (G s ]k )) ■ n jk 

5—1 


^out („, . \ _ jn* \Y\ 

“jk-y-jk ^3\^3k/J j 5 


(2.7) 
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where F = (/,#). If the fluxes are integrated with a midpoint quadrature (as suggested in [14]) 
and we use the notation u°f u jk (M jk ), uj k Uj (M jk ), Fj% := F(uJ k ), and F$ := F{uf k ), 
the semi-discrete scheme (2.7) becomes 


duj 

dt 




djk 

.in ! ^out 

L jk f U jk 


[(a; 


IPout I out 77>i 

'jk^jk d" CLj k r j k 


) ■ n jk 


■ n in n out 
a jk a jk 


(Ki - «?*)] ■ 


( 2 . 8 ) 


A basic observation that will be used below is that the semi-discrete scheme (2.8) is valid for 
any conformal grid, not necessarily triangular. All that one has to do is to make the suitable 
adjustments in the notations (e.g. |7}| being the area of the cell Tj regardless of the shape of 
that cell) and replace the sum over the three edges of the triangle by a sum over the Nj edges 
of each cell 7}. When approximating solutions to balance laws of the form (2.1) the scheme has 
one additional term due to the source term 5(tt, x, y, t), i.e. 


du : 

dt 


Ni 


w= 1 v 

IT ' 


h 


jk 


k~l 


a .. 

jk 


-.out t v J 


\(d? k Fj£ + aZ?F%.) 


M jk 


jk 1 jk/ 


a° ^ 


+SAt). (2.9) 


Here 



( 2 . 10 ) 


is a discretization of the cell-average of the source term that should be obtained with an ap- 
propriate quadrature. It is the discretization of (2.10) that serves as the topic for the next 
section. 


3 A Well-Balanced Scheme 

In this section we present a scheme for approximating the solution of (1.1) which is balanced 
via a discretized average of the source terms. In Section 3.1 we show that the KP scheme 
(2.9) cannot be balanced using a straightforward discretization of the source terms on general 
conformal unstructured meshes. In Section 3.2 we present a new scheme based on (2.7), which 
does allow such balance on general conformal triangular meshes. 

3.1 Balancing the KP scheme 

Our goal now is to look for a discretization of the source term (2.10) such that the scheme (2.9) 
will preserve stationary steady-state solutions. Hence, we assume zero velocities, u = v = 0, and 
a constant surface, i.e. w = Const. 

Since the first equation in (1.1) is homogeneous, we start by considering the second equation. 
Similar analysis applies to the third equation. We are therefore looking for a discretization of 
the average of the source term 


so -i \ 

w-v 


_ B)B X . 
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over the cell Tj. The velocities are zero and therefore the only non-zero component of the flux 
in (1.1) is 

(3.2) 

This means that in order for the scheme (2.9) to preserve stationary steady-state solutions, the 
average of the source term (3.1) over the cell T } has to be discretized such that for / given in 

(3.2), 


i Nj h 

0 = 1 y n jk 

m \ £;<>?> + a 


X jk 




(3.3) 


Here njk }X is the component of the normal Uj k in the ^-directi on. 

The eigenvalues of the Jacobian of the system (1.1) are u±\/w — B and u, which are ±yjw — B 
and 0 in the case of zero velocities. If we assume that w is constant and that the point values of 
B are known, the one-sided velocities satisfy af k = Under these assumptions the condition 

(3.3) can be rewritten as 


A- 


1^1 h 


hjk 


y 'out I fii 

jk Jj , 


jk 


W'jkiZ 


Nj 


_1 

MU 


~ [(w - B jk f rijk,x] , 


(3.4) 


where B jk — B (M jk ). 

We now' assume a discretization of the cell-average of the source (3.1) of the form 


*3 

S j = ~Y J 9 j k(.w-B jk )T>, (3.5) 

fc= 1 

w'here D « B x , and gj k are yet to be determined. To simplify the notations we denote m k = 
hj k rij^x- It is easy to check that m k — 0 and hence in order for the representation (3.5) to 

be consistent with (3.4) we must have 


1 Etoi miBjtlBjk - 2w) 

2 W 


Since D in (3.6) should not be a function of w we are seeking constants a k such that 


N 


Y^rn k B k {B k - 2m) = 


k = 1 


N 


Y^9k(w - B k 


.k = 1 


N 


^ ^ &kB k 


fc= 1 


(3.7) 


where "for simplicity we omit the obvious ' j -dependence from ail the notations. Eq. (3.7) can be 
rew r ritten as 


' N N N 

—2 E mkB k = E g k E a k B k) 

k=l k = 1 k = 1 

(3.8) 

N N N 

E mkB l = - E Yl akBk - 

V. — 1 fc=l /™1 
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Figure 3.1: An admissible triangular mesh 

The coefficients of the powers of in (3.8) produce the following system of equations: 

-2 m 1 = at EL 9k, i = 1, • • • , N, 

< rrti — — digi , i = 1, . . . , N, (3*9) 

, di9j = -a,*#, i, j = 1, . . . , W 

Finally, from (3.9) we have 

= (3.10) 

Z A:— 1 

Eq. (3.10) can generally hold only for N = 2. which means that one can expect to be able to 
balance the scheme (2.9) for stationary steady- state solutions only in cells that have two edges 
that contribute to the flux in the ^-direction. This obviously excludes most meshes. There are 
two cases of special interest: 

1. When the mesh is composed of triangles with one side that is parallel to the T-axis (see 
Fig. 3.1), each cell has only two edges that contribute to the flux in the z-direction. In this 
case the system (3.9) can be solved. This result will enable us to introduce in Section 3.2 a 
modification of the scheme (2.9) that can be balanced on general meshes. Due to symmetry 
considerations, such a mesh will not satisfy the balance conditions in the y-direction (that 
come from the third equation in (1.1)) unless all triangles are right triangles that are aligned 
with both coordinate axes. 

2. The scheme can be balanced in both directions in the very special case of Cartesian grids. 
This corresponds to the case previously solved in [12]. The results of [12] when viewed 
from the point of view of the system (3.9) amount to the equality 

B l (B 1 - 2 w) - B 2 (B 2 - 2 w) = (2 w - B } - B 2 )(B 2 - Si). 

Remark . If we assume a more general discretization of the cell average of the cell-average of 
(3.1) of the form 

N N 

Sj = - E 9k(p - B k) E a kk'B k ', 

k—l 1 


(3.11) 
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then the system (3.9) is replaced by (denoting m ft = mk/(2\Tj\)) 


2?7Tfc — y ( i-/— i Qk'Q'k'k', k 1, .... TV , 

< Tflfc = - QkCLkk , k = 1, • • • , N, 

, gk&kk' = — gk'O-k'k, k^k', k,k' = 1 , . . . , N. 


(3.12) 


Eq. (3.12) can be solved in the general case (unlike Eq. (3.9)). For example, in the case of a 
triangular grid (TV = 3) one possible solution (assuming mi ^ 0 and m3 + 2 m 2 7^ 0) is 




m3 -f 2m2 
2 raj 
2ra 2 

2ra2+m3 

1 


ra,3 \ 
2m 1 \ 

m3 I 
2m 2 +ra3 I 

-2 / 


(3.13) 


If mi = 0 the expression (3.13) can be replaced, e.g., by 


9 = 


0 . \ 

+ iH > 

a — j 

p 

m3 + 2mo 
mi 

2 m 2 

i \ 

m 3 _ | 

2m2*f ra-3 

2ra 2 +raa j 

V f J 


\i 

l 

-2 J 


(3.14) 


If m3 + 2m 2 = 0 a permutation of the numbering of the sides of the triangle will provide a 
solution. Other solutions of (3.12) exist. 

While formally being able to balance the scheme with the expressions (3.13) and (3.14), it 
remains unclear in what sense YLe a kk’Bk> approximates the derivative B x . In other words, it is 
not obvious that the consistency of this discretization can be established. We therefore consider 
this approach unsatisfactory. 


3.2 A new well-balanced scheme 

We focus on conformal triangular grids. Our ideas can be easily extended to other conformal 
unstructured meshes. In order to create a balanced scheme for (1.1), we propose to decompose 
every triangle Tj into several triangles as explained below. This requires a different decomposition 
for the two flux components. For the first component, /, we decompose Tj into two triangles, 
each of which has an edge parallel to the r-axis. For the second component of the flux we split 
Tj into two other triangles, each of which has an edge parallel to the y-axis. 

We denote the vertices of Tj as V)*,, k = 1,2, 3, and define the edges of Tj as E }i = V 0 i — Vq, 
Ej 2 = Vj 3 — Vj 2 and Ej 3 = Vj 1 — V) 3 . The midpoint of the A;-th edge of Tj is Mjk- I {a, b ) denotes 
the closed interval with endpoints a and b. We also use Vjk, x (or Vjk,y) 4° denote the x (or y ) 
component of Vjk- To simplify the treatment, "we - classify the” triangles into two categories as 
portrayed in Fig. 3.2: 

• Type 1: one vertex Vj\ bisects the opposite edge in the y direction: the y-component 
V 3 i. y £ I (' Vj 2 . y , Vj 3 , y ) and a different vertex V) 2 bisects the opposite edge in the x direction: 
Vj2,z £ I O'jl.ar, Vj3,x)~ (see Fig. 3.2 (left)). 


Type 2: The same vertex Vj\ bisects the opposite edge in both the x and y direction: 


t/.. c. T / !/.„ t/.„ 1 

r ]L,y 


3 nr? T/., 

— - ' Ji.x 


T (V, n _ 

■* \ ■ j ^ 


lAo J 

• JO.X / • 


(see Fig- 


3.2 (right)). 


9 
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Figure 3.2: The two types of triangles Left: Type 1. Right: Type 2. 



Figure 3.3: The triangle decomposition for the /-flux. 

Note that these definitions specify our vertex numbering convention. The decomposition for 
the ^-flux will depend on the type of triangle. 


3.2.1 The Decomposition in the ^-Direction 

For the flux in the x-direction, we decompose 7} into two triangles Tf' 1 and Tf' 2 (see Fig. 3.3). 
Define Vff 1 to be the intersection of edge Ej 2 with the line y = V 3 i, y - The vertices of Tf' 1 
are Vff 1 := V)i, V# 1 := Vj 2 , and V**’ 1 . The vertices of Tf' 2 are Vff 2 := Vff , Vff 2 := V j3 , 
and Vff := V}i- For both triangles Tf' 1 , the lengths of the sides are hff, the corresponding 
midpoints are Mff, the interior and exterior speeds of propagation are a 


are nff = (• 

r.l , ,i.2 


jk 


jk 


and the normals 


X.Z 

k.x 5 


hj2 + = hj2: Tl'i — TLjl.x', n 


£ ). These definitions imply Mff = Mji, Mff = Mj 3 . tiff = hji, hff = hj 3 , 

= Hj 2 .x: nfi = nj 3 ,x, and nf'j x = nff 2 x = 0. For the 

•-.X,l ; in ,out ^in.out.X.2 = a in.out_ ^ a >n,out,X.2 __ Q in.o U t_ 


X. 

jk 

x.l 


x,l 


— n 


x.2 

jl 


, r ^ • i • i in.out.X.l in. out in. out 

velocities we have a sl ' = . a„* 9 = a 


ji 


j2 


a 


"fi 


7 2 


velocity 

Finally, the external reconstructions are denoted a-s 


does not plaj" any role in the following computation since n x % x = = 0. 


03 


x.l 


x.l x.2- 

7/. =7/ .; 


7/. -o prid 7/ x F = 7/.-o 
“U** 
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On each triangle T f' t we denote the cell-average of the component of the flux in the ^-direction 
by . It is given by 


fe i=i ' 2 - 


1 2 h x ’ l n x>i 

- L_ v -Ma e 

3 \rpXd\ / j in ,X,i , out 

l 1 3 \k=i a jk a 3 k 

^D.avpr^o’p of thp r-rnrnr>nnent of the flnx over the entire cell T f is taken to be the weighted 


average 


:= 


+ \3 


fX,2 I 


101 




\Tj\ J 


Since 0 'fc n jfc = 0 we have h ?i 


x,l x,l 


n n 


h % and hX j2 nX jtx = Hence, $,• can 


a:, 2 x.2 


'j 1 ,x 


be rewritten as 


<T) — hjiUji'X 

j 101 


a fi fjT + ^r/j! <*£/ (ttj2 (Cl)) + o° 27 (^- (Cl)) 


a?, + af* 


101 


a'ji + a jT a j 2 T u j 2 

afs fjz + afffjz ^ &/ («j2 (Q 2 )) + affi/ (u,- (Q- 2 )) 


( 3 . 15 ) 


where Pj\ := M x 2 and Pj 2 := Mjfl. 


^ + ^3 


x.2 


«) n 2 + <$* 


Remark. We would like to note that the flux term (3.15) can be derived directly from (2.7) by 
changing the quadrature points on E j2 to Pi and P 2 . 

We now verify that <L,- given by (3.15) approximates —f x to second order. 

Lemma 3.1 A ssurrie a smooth reconstruction u n {x,y) in ( 2 - 4 )- Assume that the flux is linear, 
i.e., f(x, y ) = ax + by + c with constant a and b. Then = —a. 

Before proving Lemma 3.1 we consider the following geometrical lemma. 

Lemma 3.2 With M and P defined as above (suppressing j), 

,, P _ m , p m 

M 1,2 — ri ;X = —7 , M3 ,x - p2,x — T • 

^2^2.x ^2^2, x 

Proof. Mi = \ (Vi + V 2 ) , and for some 5, Pi — V2 + s (V3 — V2). We require that P i>2/ = Mi >y , 
which implies that 

1 hi y “ W.y 

5 " 2 V Sty -V 2i y 


Therefore 


Mi ?x ~ Pi : ~ = — (V 1;X + Vo. x ) — if 


1 i 1.1/ ’ // 2 : y ZT/ T7 

_ T/ r (d'3.x * 2,i 


O T7 

* v '^,y 


r ‘AV 


1 Pi x E 2 

0 ‘ 
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When the orientation of the vertices of 7} is clockwise, E x x E 2 = —2 |T| and 

n -2 — (—E 2: y, E 2 ,x) while when the orientation is counter-clockwise E x x E 2 = 2 |3T| and 

n 2 = ^(E 2 , y ,-E 2 ' X ), so 


— Pi x = 


m 

h 2 n 2 .i 


Similar arguments bold for M z ~ — P 2 _. 

Proof, (of Lemma 3.1) Since the reconstruction (2.4) was assumed to be smooth, we have 
uj k = Ujf and therefore a'f k = In this case 






| T | w t x- 1 -tyj j^oj 

where we use the notation /(P ) = f(uj{P)). For the linear flux /(a:, y) = ax + by + c we have 
— — [ail7i iX + bM X y — (aP^x + 6Pi y )] l a M 3 , x + bM 3 _ y — (aP 2x + bP 2 , y )} 


\T ! 

_ hirai. s 

^T 


m 


U Pl.x) 

\T\ 


a - 


Pl n 2,3. 


hn^x 

m 

h 3 n 3 


\T\ 


a i^M 3x P 2 .x ) 

m 

h 2 n 2 . x 


a - 




The second equality holds since M x<y = P Xt y and M 3tV = P 2 . y , while the third-equality is due to 
Lemma 3.2 and ]Pfc=i hji-Jijk = 0. 


3.2.2 The Decomposition in the y-Direction 

Type 1 triangles. We decompose 7} into two triangles Tf’ 1 and Tf’ 2 (see Fig. 3.4 (left)). The 
intersection of the edge Ej 3 with the line x = Vj 2 . x is denoted by Vjf 1 . The vertices of Tj' x are 
Vf{\ Vf 2 } := V jU and Vft 1 := V j2 . The vertices of 7f 2 are vertices Vff := V j2 , Vff := V j3 , 
and 1/// := V/f 1 . As before, we have M]f = M jU Mff = M j2 , hff = h jU hf = h j2 , 
h ji + h j 2 = h j 3 - The normals n v j2 = n jU nff = nf 2 = n j3 , nf 2 = n j2 , nff v = nff y = 0, 


L'i 


L'i 


in, out. y. 2 

J, j2 > a j2 


and the velocities a^T 1 *’ 1 = aj 3 ’ out , a^ out>y ’ 1 = aT"*, a£ 0ttt ’ v ’ 2 = 0 %"*, c&T*'* = The 


j3 


velocity does not appear in the result because nff y = nff y = 0. Finally, the external 

reconstructions are uff = Uj 3 , uff = u jU uf 2 = u j2 , and uff = % 3 - 

We denote the cell-averages of the component of the flux in the y-direction on each triangle 
T?’ 1 by r*. It is given by 


1 2 h y ,i y.i 

pi , = 1 XT n jk ri jk..y 

i ' | jV’* I 


k=l a jk 


-.vs 


Out. 

°jk 


3? Krt « (y¥)) + a jT yA 9 («, «))] . 


• = 1.2. 


The cell-average of the y-component of the flux over the entire cell 7} is then given by the 


weighted average 


\T y '- T?'* 

r f 3 r i i ! 3 i t" 2 

1 3 ■— \rr i 1 7 1 - 1 


Mj! 


I.qn I 7 * 
[TW 
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Figure 3.4: The triangle decomposition for the 5 -flux. Left: Type 1. Right: Type 2 . 


Clearly = -h)?n v £ y , and h y ?n y ? y = -h$r$ y . Therefore Tj can be rewritten as 

hjin^.y f + (LjTgfi a%9 (u j3 (Pjs)) + a$g (vg (j^s))l , 


aft + a° u f 


“ft + “ft* 


h j2 n j2 , y T aft<$‘ + a°$g% a%g (u jZ (P J4 )) + a $9 fa (Pj*)) 

‘ 12} I l “ft + “ft* “ft + “ft* 


“ft + “ft* 


where P,a := Mf ,’ 1 and P ,- 4 := 



Type 2 triangles. This case corresponds to Fig. 3.4 (right). Analogous computations to those 
for Type 1 triangles provide the cell-average of the y-component of the flux over the cell Tj, 
which this time is given by 


hji n jhy 

T oftc^f + 

af 2 g (uj 2 (Fft)) + “ft‘P (uj (Pjs)Y 

n 

L “ft + “ft* 

“ft + “ft* 

h , j3^j3 > y j 

\a% 9 °f + a^fg% 

“ftp («*2 (Pj*)) + a jfg («i ( p ji)Y 

n 1 

. “ft + “ft 1 

“ft + “ft* J 


where Pj Z and Pj 4 are given in Fig. 3.4 (right). 


3.2.3 The New Method (for conservation laws) 

We would now like to combine all the different ingredients that we developed in the previous 
section into one scheme. The scheme that we write here is still a scheme for approximating 
solutions of the conservation law (2.2) without the source term. Based on our preliminary 
analysis of Section 3.1 we know that we will be able to find an admissible discretization of the 
source terms that will result with a well-balanced scheme. We will treat the source terms in the 
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We write two versions of the scheme based on the type of triangle. For type 1 triangles, the 
discretization of the component of the flux in the x-direction, cfv is given by (3.15) while the 
discretization of the component of the flux in the y-direction, I),, is given by (3.16). In this case 
the scheme takes the form 


hjlTljljc 

‘ — in /out -out /in 

12} 1 

ft'j3'R'j3,x 

„m 1 -out 

U jl ^ a j 1 
a jzfj3 + «j3 /j3 

|2}l 1 

a; 3 + a# 

hjlTljly 

faftpjr + 

12} 1 

. a ) n i + S°i‘ 

hj2'H'j2,y 


12} | 

-in -out 

a j2 ' a j2 

1 3 
— V/? 

-in -out 

Lout 

!r-i ^ h 

^ j fc = 1 

— in 1 —out L 

T 


a j2 + a f2 


aj 2 + a°f 


x jZ ' 3 


n in _t_ n onz 
a j 3 < a jS 


(3.18) 


Here P, = P 2 = V 3 + Ifjft. P, = V, + ifjis, and P 4 = V 3 - | JjP 3 . 

For type 2 triangles we replace the discretization of Fj by the one given in (3.17) ending with 

duj _ hjiUji^ [ a) n i//r + aff/j? (Uj 2 (Pji)) + afff ( Uj (Pji))] 


n in -L n out 

a j 1 ' a jl 


af 2 + a°-? 


VW + aff /jg _ (%2 (P j2 )) + aff / frj (gg)) ' 

|2} | L a ; n 3 + S-3 «i2 + a?2 J 


<2 + a?2 


hjin jUv o^9]f + affgft Kg (2}s)) + aff g (uj (F}- 3 )) 


1 2} I L 


a“i + a j? 


a 1 /: 2 + a j 2 


(3.19) 


h^nj^y °^9jf + a$gf 3 _ af 2 9 (^2 (2}a)) + a$g (uj (P j4 )) 

|2}l l + a) n 2 + 




n lu _U /~i oux 

a i3 T a j3 


- «£] • 


In this case P,=V 2 - E 2 , P 2 = V 3 + ^E 2 , P 3 = ^2 - ^^ 2 , and P A = V 3 + ^£ 2 . 


Remark. A simple case of interest is that of coordinate-aligned right triangles. Such triangles 
can be considered to be of type 1 with n j2 . y = rijz. x = 0 and Pji = Mj 2 and Pi 3 = Mjz, which 
means that the method becomes (taking into account that £} fc=1 hjk n jk = 0) 


1 f— in /out I -out fin -in f out 1 -out fin 

1 r fo /ji A n . + Vi2 + fji . 

2 } 1 L + a jT a ?2 + a j2 

1 \a^g°f + affgfr t , a;" 3 p/ 3 ut + ^ 

T-l n in , -out 1 -in ^ -out A j'3 n j3..» 

^1 L a 7l^ a il a j3 ‘ a j3 




a, n o a 1 


fojzTljZ'. 


■J3 1 jZ 


l -Z h ^ 


in j_ n out 
jk ' ^f/c 


o, out — 'll 111 

Pi* u J k J 


(3.20) 





14 


S. Bryson and D. Levy 


As expected, in this case the method (3.20) coincides with the KP method (2.8). 


3.2.4 Adding the Source Term 


We return to the shallow water equations (1.1). Our goal is now to find an admissible discretiza- 
tion of the cell-average of the source terms that will preserve stationary steady states. Such a 
term will be added to the RHS of (3.18) or (3.19) depending on the type of the triangle. In 
the following we assume Type 1 triangles. Similar analysis holds for Type 2 triangles. We also 
note that the first equation in the system (1.1) is homogeneous. This means that we only have 
to consider the remaining two equations when dealing with the source terms. We recall from 
Section 3.1 that in stationary steady states all the velocities are equal, aj k = a°^ = a ]k - Our 
new method (3.18) for the last two equations of (1.1) then becomes 


0 = 


U'j 1 Tlj 2 

"m 

*^3 3^1 3 .x 

‘ Wi\~ 

lij\ Tlj J y 
‘ 2 12} | 
foj2'Uj2,y 
' 2 12} [ 


( (w (Mji) — B (Mji)) 2 

V 0 

(w (Mjz) - B (Mj 3 )) 2 
0 

0 

( w (Mji) — B (Mji)) 2 
0 


(wiP^-BiPj i)) 2 y 

V 0 J. 

{w{P j2 )-B{P j2 )f 
0 

0 

( w (Pjs) ~ B (Pjz )) 2 

0 

(w(P 34 )~B(P j4 )) 2 


(3.21) 


+ ( S 


The last term in (3.21) represents the average of the source, i.e M Sj = avg(— (w — B)B X ) and 
S'? = avg(— {w — B)B y ). where both averages are taken over Tj. 

We use (3.21) to determine the admissible discretizations of the cell-averages of the source 
terms. For constant w the source terms (given by (3.21)) can be rewritten as 


% 


S ? = - 


~ B Mjl + w Pj , - B Pjl )[B Mjl - Bp n ) 
•~ 2j ~ Bmj3 + W P]2 - Bp j2 ][B Mj3 - B Pj2 ] , 

ol ' v \wM n - B Mj i + wp j3 — B Pj3 ][B Mjl — B Pj3 ] 
[wMi 2 ~ B m 2 + wp 4 — B P }[B m, i 2 ~ B Pji }. 


2im 

hj2Uj 


ji-.y 


(3.22) 


Here, we use the notation Wm 1 := w(Mji), etc. We take the expressions in (3.22) as the dis- 
cretization of the source even when w is not constant. Similar expressions can be easily written 
for Type 2 triangles. 


Lemma 3.3 The source discretizations Sj given by (3.22) are consistent approximations of the 
source terms in (1.1), and lead to detailed balance in the stationary steady-state case (3.21). 


Proof. We show that S 7 - 


TT -*™ -to-i T owi i 

■L 1 


~ avg(— (w — B)B X ). To simplify the notations we suppress the index 
mow tlW fl-L — P. - = — — .Since AT .. — P -\ we have 

w * -L.m /12T12.X ' *** 
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fi* 1 zpf 1 — B x + O ( \Mi — Pi| 2 ) at the midpoint between Mi and P\. Hence the first part of Sj 
in (3.22) becomes 


2 h 2 ri2,x 


— Bmi + w Pi 


-B Pl ) 


Bju \ — Bp x 
M hx - P ltX 


hyn^x 

2/l-27T'2.x 


(wmi — 4 - wp 1 — Bp } ) 


Applying a similar argument to the second term of Sj in (3.22) gives 


S} 


1 f h\TL i x . 

nss ( ” Mi 


Bmi + inpj — ) + 




(u)m 3 - B m 3 -b Wp 2 - Bp 2 )} B x . (3.23) 


Clearly, the coefficient of B x in (3.23) is a discretization of a weighted average of — (w — B ). 
For example, when — i?is constant we have 




h x rii~ ~t~ 
h 2 n 2iX 


(w - B) B x = - (w - B ) B x . 


Similar arguments hold for Sj. 


4 Numerical Examples 

The scheme developed in Section 3.2 did not assume any particular reconstruction. There are 
several different second-order reconstructions on triangular meshes that are being used in the 
LfPT-Qf/nrp (qpp ML! and the references therein). We briefly describe the one we used in our 
simulations. 

The starting point is the limited least-squares estimate of the gradients as done in [2]. The 
first step is to compute a least-squares estimate of the gradient of a field / on the triangle Tj , 
Vjf. We then limit the gradient T > i / component by component as 

= MM (vy tjj, V i2 /, Vjaf) , 

where Vj k f is the least-squares gradient estimate on Tjk and MM stands for the usual MinMod 
limiter 

{ min, {xj} , if Xj > 0. Vj, 
max, {a:, } , if x, < 0, Vj, 

0, otherwise. 

We use the gradients D, to construct a piecewise linear reconstruction for the point-values of 
each triangle edge E lk as 

Uj(x) = u.j -i- MM (D ju • (x — Sj) , BjkU ■ (x — x,)) . 

Here T>j k u is the limited gradient estimate on Tjk, x, is the center of Tj and x 6 Ej k ■ We find 
that this double use of the MinMod limiter minimizes spurious oscillations while preserving the 
second-order accuracy of the reconstruction. 
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Figure 4.1: The two types of triangle meshes, 
triangulation. 



We use an adaptive time step given by 
At — 0.9 min 7-^-7, 

i UU 


where ry- is the radius of 7} and jA^| is the norm of the largest eigenvalues of the Jacobian of 
the flux on Tj . 

We use two types of meshes in most of our simulations. The first is a uniform triangulation 
of a Cartesian mesh with N x N nodes and constant spacings Ax and Ay. which divides each 
Cartesian cell into four triangles. To test our method on more general meshes, our second type 
of mesh is generated from a uniform triangular mesh with a perturbation of the coordinates of 
the interior vertices. Examples of both meshes are shown in Fig. 4.1. The only exception is 
Exaniple 6, which uses a uniform triangulation of a warped Cartesian mesh. 


Example 1; accuracy tests for a 2D linear advection equation 

We test the accuracy of our method (3.18)~*(3.19) wdien applied to the two-dimensional linear 
advection problem 


u t + u x + u y = 0, (x, y) E [0, l] 2 , 

u (x, y. t = 0) = sin(27rx) + cos(27ry), 


(4.1) 


with periodic boundary conditions. The relative I 1 -error (i.e. the U-error divided by the L l ~ 
norm of the exact solution) at T = 1 is shown in Table 4.1 for both uniform and perturbed 


01 ICbllX U.1C101VJ11O 


fVi 


fi o Drvmr 


.n-rrlnr a r-r'i 


irarv of rmr method. 
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Uniform Triangulation 

Perturbed Triangulation 

N 

relative L 1 -error 

L l -order 

relative L l -error 

L 1 -order 

10 

0.123 

- 

0.129 

- 

20 

0.029 

2.09 

0.034 

1.94 

40 

6.52 x 10' 3 

2.15 

8.68 x 10~ 3 

1.95 

80 

1.73 x 10" 3 

1.91 

3.02 x 10“ 3 

1.52 


Table 4.1: ZT-error and convergence rates for the advection equation (4.1) at T = 1 on a uniform 
and a perturbed triangulation of an N x N Cartesian mesh. 


Example 2: accuracy tests for a 2D Burgers equation 

We continue by checking the accuracy of our method (3.18)-(3.19) on non-linear problems by 
applying it to the two-dimensional Burgers equation 


j u t + \ (« 2 ) x + I ( u 2 )y = 0, [-27 r, 2tt} 2 , 

1 u(x,y,t = 0) = sin (^) , 


(4.2) 


with periodic boundary conditions. The L 1 -error of our method at T = 0.5 is shown in Table 4.2 
for uniform and perturbed triangulations. The computed cell-averages after singularity formation 
at T = 1.5 are shown in Fig. 4.2. Note the sharp shocks that are captured by our method. 



Uniform Triangulation 

Perturbed Triangulation 

N 

relative L 1 -error 

L 1 -order 

relative L 1 -error 

L 1 -order 

10 

0.064 

- 

0.064 

- 

20 

o.oi4 

2.20 

0.014 

2.20 

40 

4.68 x lO" 3 

1.58 

4.65 x 10 -3 

1.59 

80 

1.37 x 10~ 3 

1.77 

1.36 x lO" 3 

•1.77 


Table 4.2: L-Cerror and convergence rates for Burgers equation (4.2) at T = 0.5 on a uniform 
and a perturbed triangulation of an N x N Cartesian mesh. 


Example 3: accuracy tests for 2D systems 


To test our method (3.18)-(3.19) on a system of equations we apply it to the Sod problem for 
the Euler equations 


P 1 


pu ' 
pv 

E 1 

4- 


/ t 



pu 

\ 


( 

pv 

\ 


pu 2 +p 




puv 



puv 




pv 2 + p 


\ 

u (E -j- p) 

i 

/ 

•V* 

i 

\ 

w 1 rj 

J 


= 0 . 


(4.3) 
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Side view 


-no 


1 o 


Figure 4.2: Burgers equation at T = 1.5 on a uniform triangulation of a 20 x 20 Cartesian mesh. 
Left: oblique view. Right: side view. 


Here p is the density, ( u , v) is the velocity, and E is the energy. The equation of state for the 
pressure is p = (7 — 1) \E — f ( u 2 ~f u 2 )] . We set 7 = 1.4 and take the initial conditions (constant 
in the y-direction) 


f (1.0, 1.0, 0.0, 0.0), x < 0.5, 
\ (0.1,0.125,0.0,0.0), x > 0.5. 


(4.4) 


Figure 4.3 shows the computed cell-averages at T = 0.16 using our method from Section 3.2 
projected onto the y = 0 plane. A reference one-dimensional solution is also shown. The two- 
dimensional problem uses a uniform triangulation based on a 300 x 30 Cartesian mesh on the 
domain [0, 1] x [0, 0.1], The one-dimensional reference solution is computed using the second-order 
central method of [13] with 3000 points in the domain [0, 1]. 


Example 4: a balance test 

In this example we demonstrate the well- balanced property of our scheme (3. IS)— (3. 19) with the 
source discretization (3.22). 

As a simple test of balance, we consider the shallow water problem with initial condi- 
tions that represent a stationary steady-state. We choose w(x, yft — 0) — 2, u (x, y, t = 0) — 
v (x. y, t — 0) =0 and a bottom topography given by B (x, y) = sin(27rx) + cos(27ry) on the do- 
main [0, l] 2 . We assume periodic boundary conditions. This initial value problem has the trivial 
stationary steady state solution w (x, y, t) = 2, u ( x , y, t) — v (x, y,t) — 0 for all t . 

The relative I^-errors at T — 1 on both uniform and perturbed triangulations based on a 
10 x 10 Cartesian mesh are given in Table 4.3. This table also shows the result found using the 
KP method [14] with the same source discretization (3.22) which is not balanced for this scheme, 
yy 0 500 that our method balances to machine accuracy. 
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O 2D triangular simulation 
1 D reference simulation 



Figure 4.3: The pressure (left) and u - velocity (right) fields of the Sod problem at T = 0.16 on a 
uniform triangulation. The circles show a projection onto tne y = o piane oi cue solution oi tne 
2D problem on a triangular mesh. The line shows a ID reference solution computed with the 
second-order central- upwind method of [13]. 


Method 

Uniform Triangulation 

Perturbed Triangulation 

(3.18)-(3.19)-(3.22) 

3.5 x lO"- 17 

6.9 x 1CT 19 

[14]..+. (3.22).- 

.. - ... 3.0 x 1(H-.- - 

3.5 x 10~ 3 - - 


Table 4.3: The relative LDerror in the balance for (3.18)-(3.19)-(3.22) and for the KP scheme 
with the source discretization (3.22) 
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Example 5: propagating waves with a bottom topography 

We apply our method to a test problem from [16] of a small perturbation of a steady state 
problem on the domain [0, 2] x [0, 1] with periodic boundary conditions in y-direction. The 
bottom topography is the elliptical Gaussian mound given by 

B (x, y) = 0.8 exp (—5 (x — 0.9) 2 — 50 (y — 0.5) 2 ) , 

and the initial conditions are 

, f (1.01,0.0,0.0), if 0.05 <x <0.15, 

w,u,v | (1.0, 0.0, 0.0) , otherwise. 

Fig. 4.4 shows the result of our method at various times on a uniform triangulation of both a 
200 x 100 and 400 x 200 Cartesian mesh. These results are in good agreement with other methods 
on Cartesian grids (see [12, 16]). 

Example 6: converging-diverging channel with bottom topography 

Our final example is that of a converging-diverging channel with critical flow adapted from [10]. 
The channel is defined on the domain [0, 3] x [—0.5, 0.5] with a half-cosine constriction centered 
at x = 1.5. The mesh for this example is shown in Fig. 4.5 (a). It is a uniform triangularization 
of the warped Cartesian mesh defined by the mapping ( x , y) — > ( x , (1 — 0.2 cos 2 (n ( x — 1.5))) y) 
when \x — 1.5| < 0.5. The initial data is w = 1, u = v = 0. The y-boundaries are reflecting. The 
left x-boundary is an inflow boundary with u — 5.0 and the right x-boundary is a zeroth-order 
outflow boundary. We run the simulations on a 90 x 30 mesh until T = 7 after the steady state 
is achieved. 

We first present this example with a flat bottom, with contours of w shown in Fig. 4.5 (b). 
Fig. 4.5 (c) shows the same channel at the same time with bottom topography 

B (x, y) = 0.8 (exp (-10 (x - 1.9) 2 - 50 (y - 0.2) 2 ) + exp (-20 (x - 2.2) 2 - 50 (y + 0.2) 2 )) . 

This topography is shown in Fig. 4.5 (d) and represents two elliptical Gaussian mounds just 
down-flow from the constriction. 
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